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Abstract 

We analyze the critical properties of the three-dimensional Ising model with linear parallel ex¬ 
tended defects. Such a form of disorder produces two distinct correlation lengths, a parallel cor¬ 
relation length £|| in the direction along defects, and a perpendicular correlation length £j_ in the 
direction perpendicular to the lines. Both and £_i_ diverge algebraically in the vicinity of the crit¬ 
ical point, but the corresponding critical exponents v \i and is± take different values. This property 
is specific for anisotropic scaling and the ratio v\\/v±_ defines the anisotropy exponent 8. Estimates 
of quantitative characteristics of the critical behaviour for such systems were only obtained up 
to now within the renormalization group approach. We report a study of the anisotropic scaling 
in this system via Monte Carlo simulation of the three-dimensional system with Ising spins and 
non-magnetic impurities arranged into randomly distributed parallel lines. Several independent 
estimates for the anisotropy exponent 8 of the system are obtained, as well as an estimate of the 
susceptibility exponent 7 . Our results corroborate the renormalization group predictions obtained 
earlier. 

PACS numbers: 05.10.Ln, 64.60.F-, 75.10.Hk 
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I. INTRODUCTION 


The effect of structural disorder on criticality remains one of the most attractive issues 
in condensed matter physics.- Realistic systems always contain some imperfections of their 
structure. Thus the question of how disorder does influence the critical properties of a 
system deserves considerable interest. Obviously the modifications introduced depend on 
the amount of disorder as well as on the type of disorder. Quenched disorder is usually 
studied in the form of dilution (random site^ or random bonc£ systems), or as a random 
field,- random connectivity^ or anisotropy.-^ 

In our study we will focus on the effect of weak disorder in the form of random sites. In 
the simplest case, such disorder may be considered as uncorrelated, randomly distributed 
point-like defects. The relevance of point-like disorder for critical behaviour is predicted by 
the Harris criterion:- it changes the universality class of d -dimensional system if the heat 
capacity critical exponent a pure of the corresponding homogenous (pure) system is positive. 
Since in d — 3 the pure Ising model has a pure > 0, weak point-like disorder there leads to 


a new critical behaviour. Results of ana 

f. 2 


ytical, numerical, and experimental studies of this 

celebrated system are reviewed in Ref. 

Many real systems contain more complex forms of disorder, for instance dislocations, 
disordered layers, grain boundaries, cavities or other extended defects. To take this into 
account Weinrib and Halperin have proposed a model,- in which defects are correlated with 
a correlation function decaying with distance x according to a power law: g(x ) ~ x~ a . 
There is possible interpretation for integer value of a: the case a = d — 1 (a = d — 2) 
describes lines (planes) of random orientation, while a = d corresponds to above-mentioned 
uncorrelated defects. The critical properties of such systems have been extensively studied 
by the renormalization group (RG) approach^— as well as through Monte Carlo (MC) 
simulations^ - —. 


Another possible implementation of extended defects was proposed by Dorogovtsev^ 
within the model of a d -dimensional spin system with quenched random impurities that are 
strongly correlated in dimensions and randomly distributed over the remaining d — 
dimensions. On the contrary of the model of Weinrib and Halperin which is isotropic, - 
Dorogovtsev’s model describes a system which is expected to behave differently along the 
directions “parallel” to the ^-dimensional impurity and in the “perpendicular” hyper-planes. 
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The case £d = 0 corresponds to point-like defects, and extended parallel linear (planar) de¬ 
fects are respectively given by £d = 1(2). Generalization of £d to nonnegative real numbers 
may be interpreted as an effective fractal dimension of a complex random defect system, 
ffowever, the relation of analytically continuation to non-integer Euclidean dimension and 
a fractal dimension is not straightforward.— Critical behaviour of these systems was exten¬ 
sively studied by RG methods with the help of double e = 4 — d, e^-expansions^^ as well 
as applying resummation technique to the RG asymptotic series directly at d = 3 and fixed 

e J_9-_22 

There exist also investigations in the mean-field approximation^ as well as MG simula¬ 
tions which have some connections with Dorogovtsev’s model with Ising spins and extended 
defects with = 2.—' — In these latter studies, disorder was modeled by random bonds 
between planes of spins. To our best knowledge, the model with parallel linear extended 
defects in dimension d = 3 has not been studied numerically so far, although the existing an¬ 
alytical predictions are accurate enough to challenge MC verification. Therefore we present 
such MC study in our paper. The aim of this paper is therefore precisely to report such an 
extended numerical study. 

The paper is organized as follows. In the next section we briefly present the analytical 
predictions for scaling of three dimensional (3d) system with parallel extended defects as well 
as the RG estimates for the critical exponents of the Ising system with randomly distributed 
parallel linear defects. In Section IIIII we remind the essentials of anisotropic finite size 
scaling. The formulation of our model, definition of the observables as well as details of 
the simulation are listed in Section IIV1 We eventually present the results of simulations in 
Section [V] and Section ED summarizes our study. 


II. ANALYTICAL RESULTS FOR A SPIN MODEL WITH EXTENDED DE¬ 
FECTS 


The critical behaviour of the model under consideration is described by the following 
effective Hamiltonian: 


n 



M + V(x)f‘ 2 (x) + (V ± 0(x)) 2 


+« O (V||0(x)) 2 ) + ^y(0 2 (x)) 2 


( 1 ) 
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Here, </> is an m-component vector field: 0 = {0 1 • • /io and Uq are the bare mass and 

the coupling of the magnetic model, cco is the bare anisotropy constant, and V(x) represents 
the impurity potential, which is assumed to be Gaussian distributed with zero mean and 
correlator: 


V{x)V{x') — —v 0 5 d £d (x± — x'±). 


( 2 ) 


Here, overline stands for the average over the potential distribution, (-fo) is a positive 
constant proportional to both the concentration of impurities and the strength of their po¬ 
tential. The impurities are envisaged as ^-dimensional objects, each extending throughout 
the system along the coordinate directions symbolized as x\\, whereas in the remaining d — £d 
dimensions they are randomly distributed. Operators Vj_ and V|| mean differentiation in 
the coordinates x_l and X||, correspondingly. One assumes that the linear size of the defects 
is much larger than the spin-correlation length and also larger than the linear separation 
between any pair of defects. This assumption is valid for defect concentrations well below 
the percolation threshold. 

The static critical behavior—^— of this model as well as critical dynamics near 
equilibrium^ 1 ^ - — were examined by means of the RG method. A double expansion in 
both e, £d was suggested and RG functions were calculated to the first order.— These results 
were consistent with a crossover to a new universality class in the presence of extended de¬ 
fects. These calculations were extended to the second order in Ref. [l8l Here, it was argued 


that the Harris criterion is modified in the presence of extended impurities: randomness is 
relevant, if 


2 

£ d > d -, (3) 

^pure 

where z/ pure is the correlation length critical exponent of the pure system. For point-like 
defects (ed = 0) Eq. (J3J) reproduces the usual Harris criterion. Resummation technique 
applied to the RG asymptotic series^ - — lead to reliable estimates of the critical exponents 
for this model. Furthermore, different scenarios for the effective critical behaviour were 
discussed.— Within RG approach, the influence of cubic anisotropy of the order parameter— 
and the effect of replica symmetry breaking on the disorder average^ were studied. Reports 
of short time critical dynamics in systems with extended defects are also available in the 
literature.— For completeness, we also mention here several papers more, where models 
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with more complex forms of disorder, including extended defects as a particular case, were 
analyzed.— ^ 

The model described by the effective Hamiltonian (jT]) has rich scaling behavior. As it 
was already mentioned, such a system is no longer isotropic. Due to the spatial anisotropy, 
two correlation lengths exist, one perpendicular and one parallel to the extended impurities 
direction: (£j_ and £||). As the critical temperature T c is approached, their divergences are 
characterized by corresponding critical exponents u j_, v\\. 

~ \t\~ v± , £n ~ \t\~ un , (4) 

where t is the reduced distance to the critical temperature t — (T — T c )/T c . The corre¬ 
lation of the order parameter fluctuations in two different points acquires an orientational 
dependence.— Thus, the critical exponents r/± and 7711 , that characterize the behavior of 
the correlation function in the directions, perpendicular and parallel to the extended de¬ 
fects, must be distinguished. Spatial anisotropy also modifies the critical dynamics near 
equilibrium producing two dynamical exponents z\\ and z±.— On the other hand, as far 
as the interaction of all order parameter components with defects is the same, the system 
susceptibility is isotropic^ and can be expressed by the pair correlation function:— 

X(k±,k\\,t) = 


In (J5]), k\\, k± are the components of the momenta along £ d and d — Ed directions, respectively, 
7 is the magnetic susceptibility critical exponent and g is the scaling function. The following 
scaling relations hold:—:— — 


ltl-7,, | _t±_ _h_ ± i 
11 9 1 


k v± ~ 2 n ( 1 1*1 

K± 9 \ A fe q|/-_L> k V«± 


k VU ~ 2 n 1 k 

K \\ 9 


v ±/ v 


-, 1 , 


1*1 


(5) 




7 = (2 — T)±)vj_ = (2 — 7 ||)^||. ( 6 ) 

The critical exponent a of the specific heat is related to u\\ by the hyperscaling relation 
that differs from the ordinary one:— 

a = 2-(d-£ d )v ± -£ d v\\. (7) 
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All the other scaling relations are of the standard form. This implies that one should 
calculate at least three independent static exponents (e.g., v\\, 7 ) instead of two, as in 

the standard case, to find the others static exponents by scaling relations. 

We are mostly interested in the case of Ising spins rn — 1 and linear parallel extended 
defects = 1. For this case we list below the numerical estimates of critical exponents 
obtained within RG approaches. These estimates based on yfe expansions are the following:— 

v j_ = 0.67, u\\ = 0.84, 7 = 1.34, z± = 2.67, z\\ = 2.22. ( 8 ) 

More reliable estimates are obtained with help of resummation of two-loop RG functions for 
static exponents:— 

v± = 0.750, u\\ = 0.880, 7 = 1.483, (9) 

and for dynamic exponents:— 

z x = 2.418, z\\ = 2.217. (10) 

It is desirable to check these results with MC simulations. In the following section we show 
how finite size scaling may be used to extract quantitative characteristics of the critical 
behaviour from MC data obtained for finite systems. 

III. ANISOTROPIC FINITE SIZE SCALING 

Systems studied by MC simulations have limited sizes. Therefore finite size scaling 
(FSS)^iF— plays a key role to extract the critical exponents of thermodynamical functions 
from MC data. According to FSS theory, in the vicinity of critical point t —> 0 the order 
parameter M and the susceptibility x of isotropic systems scale with the linear size L as: 

M(t -A 0) = L-^ U M(L), (11) 

0 ) = L^x(L), ( 12 ) 

where A4, x are magnetization and susceptibility scaling functions. A quantity of special 
interest often used within MC simulation, namely of the fourth-order Binder’s cumulant,— 
obeys the following behavior: 

m, L) = 1 - Aft - = UiitL 1 '"), (13) 
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where ( M 2 ), (M 4 ) are second and fourth moments of the distribution of order parameters, 
C /4 is the scaling function. In isotropic models all curves of the temperature dependence 
of C /4 for different system sizes intersect in one point, indicating the location of the critical 
temperature. It is true even if one keeps different aspect ratios for the simulation box: the 
curves for different sizes at a given generalized aspect ratio all intersect in the thermodynamic 
limit at the critical temperature.— 

There exist different systems that manifest strong anisotropic scaling at criticality, i.e. 
the critical exponents for their correlation lengths differ in different directions. Among them 
we mention as examples the behaviour of the next nearest neighbor Ising model at the 
Lifshitz point ,- 0 non-equilibrium phase transition in driven diffusive systems,— dipolar in¬ 
plane Ising model,— driven Ising model with friction,— Ising model under shear,—— interface 
localization-delocalization problem.— A common feature of the above mentioned systems is 
the presence a single anisotropy axis which results in the existence of two characteristic length 
scales: L\\ along the anisotropy axis and Lj_ in the perpendicular directions. According to 
the generalization of the FSS concept for systems with two characteristic length scales,— 
the properties of the thermodynamical functions depend on the “generalized aspect ratio” 
L\\/L e ± , with an anisotropy exponent 9 = v\\fv±_. Therefore relations (fTTTh (fT 2 j) are modified: 


M(t^O) = (14) 

x(t4 0) = l/'jilln/ll). (15) 

Note that (TT4f) . (fl5|) can be equivalently represented as: 

M(«-M1) = L-/ l ^M(LJL\( e ), (16) 

X(i-X>) = Ll^xtWA"). (17) 

Anisotropic finite size scaling for Binder’s cumulant at the critical point predicts 

U A (t -► 0) = U^Lw/L^). (18) 


Therefore, one has to keep generalized aspect ratio L\\/L e j_ fixed while performing MC 
simulations to extract critical exponents. One can use the values of the exponents (151) - (fIUD to 
evaluate the anisotropy exponent 6 for the model under consideration. Using results of ^fe 
expansions (jHJ) we can estimate 9 = v\\/v± — 1.25. Taking into account that v\\/v_l = z±/z\\^ 
we can get the second estimate of the yfe expansions: 9 = z±/z\\ — 1.20. Next estimates can 


be obtained with help of the results of resummation of two-loop RG functions: 9 = v\\/i'± = 
1.173 from (jUJ) and 9 = z±/z\\ = 1.091 from (fTOlh All estimates give 9 > 1. This reflects the 
physical situation: in parallel directions, the fluctuations are stronger because they are not 
limited by defects and the correlation length in this direction diverges more sharply. 

In order to obtain an estimate of 9 directly from the MC simulations, independently of 
the RG predictions, we can make use of the relation that characterizes systems possessing 
anisotropic scaling when L± L\\ at t —> 0:—^ 

£l(£||) ~ L\(‘. (19) 

Such a relation was successfully applied for the analysis of the second order phase transitions 
in Ising models with frictions^ and under shear.— 

In the following section, we use the formulae (fTT]) -([T9 ]) to perform the analysis of the 
results of MC simulations. 

IV. DETAILS OF THE SIMULATIONS 

We consider the 3d Ising model on a simple cubic lattice of L x x L y x L z = V sites. Each 
site i = ( x,y,z ), 1 < x < L x , 1 < y < L y , 1 < z < L z , is characterized by the occupation 
number q = { 0 , 1 } with q = 1 if the site % is occupied by a spin s* = ± 1 , while q = 0 
corresponds to a non-magnetic site %. Nearest neighbour spins interact ferromagnetically 
with constant exchange J > 0. Explicitly, the Hamiltonian of the model may be written 

77 — CiCjSjSj, (20) 

where the sum (i,j) is taken over the pairs of nearest-neighbor spins. We use periodic 
boundary conditions. 

The non-magnetic sites (q = 0) in our model have to be arranged into parallel lines and 
oriented along a given axis, say the z axis. The length of these lines coincides with the 
system size in z direction, L z . Fig. Q] shows one possible configuration of such lines. 

We generate the impurity distribution for a given impurity concentration p for a lattice of 
size V in the following way. We compute the integer number n imp = int (pL x L y ) of impurity 
lines of length L z . We distribute randomly intersections of these lines with the x — y plane, 
as shown in Fig. [Q 


9 


X 


FIG. 1: (Color online) An example of the lattice of size 4x4x4 with Ising and (dark circles) 
non-magnetic impurities (light circles, red online) collected into ni mp = 4 lines. 


We compute powers of the magnetization (k = 1,2,4) 

M ‘ = ((E«)* < 21 ) 

Id 

where {i} means the lattice summation over all sites. From powers of the magnetization we 
construct the magnetic susceptibility and Binder cumulant in the following way 


x = pv((m 2 )-(\m\) 2 ), 


U, = 1 


_ ( W 


( 22 ) 

(23) 


\3(M 2 ) 2 

Here at first step we perform the thermal averaging of the power of magnetization for each 
impurity realization (denoted by angle brackets, e.g. ( M 2 )). Then we average these results 
over different impurity realizations. The averaging over impurity realizations is denoted by 
the over line. In our simulation we fix J/fc# = 1, in which case (3 = 1/T is the inverse 
temperature. 

We compute also the correlation lengths for directions parallel to impurity lines £|| and 

perpendicular to them £j_ with the help of the Fourier transform in a similar way as it 

was done for isotropic systems.— Let us introduce s(k) = e^CjSj where the wave 

j 

vector is k = (k x ,k y ,k z ) and we denote (xo) = (|s(0,0,0)| 2 ), (x?) = (\s(2tt/L x , 0,0)| 2 ), 
(Xi) = (|s(0, 0, 27 t/Z/~)| 2 ) thermal average of Fourier components for particular realization 
of impurities. Then we can compute the corresponding correlation lengths averaged over 
impurity realizations 

1 


£l= 


(Xo) 

2sin(7 t/L x )]I (xf) 


C||= 


(Xo) 

2 sin(7r /L z ) V (xf) 


- 1 . 


(24) 
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The simulation is performed with hybrid MC method.— Each of MC step consists of one 
flip of Wolff cluster followed by E/4 attempts to flip spins in accordance with Metropolis 
rule. 

In our study we deal with a concentration of impurities p = 0.2. Such concentration is 
very often used in MC simulations of Ising model with disorder,—”—^ since in this case 
concentration of magnetic sites 1 — p = 0.8 is far from the percolation threshold and from 
the pure system. Additional empirical reason to take this value of p is the following: for 
3d Ising model with uncorrelated impurities, correction-to-scaling terms were found to be 
minimal at concentration of magnetic sites O.8.— 

V. RESULTS AND ANALYSIS OF SIMULATIONS 

In this section, applying the above described formalism, we give three estimates of the 
anisotropy exponent obtained from anisotropic FSS predictions for different quantities. From 
the FSS of susceptibility we get also an estimate for 'y/v±. 

A. Computation of the correlation length 

To use relation dll?]) we need to perform calculations at the critical temperature of an 
infinite system. Let us assume that the value of the critical temperature of an infinite 
system is very close to the value of the pseudo-critical temperature of the system with 
maximal size. In this case that is the temperature of the susceptibility maximum y max . We 
perform the computation for a cubical system of size 128 3 with concentration of impurities 
p = 0.2. Disorder average may be performed in either of two possible protocols: In the first 
procedure, at each disorder realization, the temperature dependence of the susceptibility is 
obtained. Then all these curves are averaged to get a single curve of x that depends on /3. 
Then the value f3 at which y(/3) has a maximum is associated with the critical temperature 
(we refer this way as method A ). The second alternative is to average temperatures at which 
Xmax is achieved for each disorder realization (we refer this way as method B). 

The results obtained with both methods are given in Fig. [21 The two curves displayed 
by continuous lines in the figure show x(/3) behaviour for two separate samples, whereas the 
dotted curve shows the y(/3) averaged over 320 samples. Using the method A we first average 
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FIG. 2: (Color online) Maxima of magnetic succeptibility x ftmax = 0.2565 for method A and 
Anax = 0.2572 for method B for a system 128 3 and p = 0.2. 



FIG. 3: (Color online) Results of the fit of £ as a power law dependence of L z for a system of size 
128 x 128 x L z and p = 0.2 at points /S^ax = 0.2565 for method A and (3 ^ ax = 0.2572 for method B. 

the susceptibility and then find the point /3^ ax = 0.2565 - magenta diamond in Fig. El Using 
the method B we first find a maximum of susceptibility for each disorder configurations, and 
then average location of these maxima over samples = 0.2572 - black square in Fig. El 


With the value of the critical temperature at hand we can perform a study through 
relation (fUl to extract the anisotropy exponent 6. To this end, we analyze the correlation 
length £j_ of the system of size 128 x 128 x L z varying L z for both values of /3 max . Results 
of simulations for L z = 25 — 50 are given in Fig. El where we plot £ for _L and || directions 
as function of L z . 
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Then we perform the fit of data for £j_ in accordance with formulae: 


Z±(L g ) = aL\ , (25) 

using fitting parameters a, b. For method A we get a = 1.27(4), b = 0.76(1), while for method 
B we estimate a = 0.91(3), b = 0.90(1). Comparing (125]) with (fT9]l we get 9 = 1/6 with 
the following estimates: 6 « 1.31 ( method A) and 6 ^ 1.11 ( method B). Taking these values 
as an accuracy interval of 6 determining, 1.1 < 6 < 1.3 we see that they are in reasonable 
agreement with existing analytic RG estimates (see Section m- 

B. Computation of the Binder cumulant 

Another way to identify the right value of the exponent 9 is the following. We expect 
that all curves for the temperature dependence of the Binder cumulant I/ 4 for systems of 
different sizes but fixed generalized aspect ratio with proper 6 will intersect at the critical 
point. Therefore we can consider the system of size L x x L y x L z with L x = L y and keeping 
L z = L e x for various values of 6. For this system we compute the cumulant f/ 4 (/3,L x ) as a 
function of the inverse temperature /3 for various values of L x . At the proper value 6* we 
expect the intersection of graphs. We perform simulations for L x = 20,40,60,80,100, the 
averaging is performed over N- imp = 128 realization of impurities. 

In Fig. 01(a)- (d) we plot the cumulant for values 9 = 1, 1.1, 1.2, 1.3, respectively. Each 
curve corresponds to I/ 4 (/3) averaged over 128 disorder realizations. As one can see, the 
presence of disorder smears the single crossing point into a region of temperatures where 
all curves cross. The narrowest region is expected for the value of 9 which is sufficiently 
close to the real one. To analyze this situation for different values of 9, we use the following 
procedure. We split the total set of A r imp = 128 different impurity realizations onto eight 
series, and compute the average value (averaged over 128/8 = 16 realizations) of the invariant 
f/ 4 for each set k = 1,2,... ,8. Later on we use the averaging over U% for evaluation of 
numerical inaccuracy. In Fig. [51(a) we plot for comparison results for U% as a function of 
/? for two different series k = 1 (lines) and k = 2 (triangles) for 9 = 1 and various sizes 
L x = 20, 40, 60, 80, 100. In this figure we observe an important difference between graphs 
for two series due to impurity induced fiuctuations. Our aim is to study the scattering of 
intersection points. Ideally, graphs for five system sizes should intersect in ten points. In 
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FIG. 4: (Color online) Magnetization cumulant U± as a function of the inverse temperature (3 for 
L z = L e x . (a) 0 = 1.0; (b) 9 = 1.1; (c) 9 = 1.2; (d) 9 = 1.3. 


Fig. [51(b) we plot data for U^" 1 and indicate the intersection points (/?*), U^) by black circles. 
Here a pair ij of indexes labels two system sizes i — 1,2, 3, 4, 5 for L x = 20, 40, 60, 80, 100. 
The graphs for some system sizes have intersections outside of considered interval for /?, e.g., 
for L x = 40 and for L x = 80, 100. In this case we select as “intersection point” the point 
between left side or right side ends of these two lines selecting the side with the smaller 
distance between the endpoints. 

Another possible situation with multiple intersection points may happen due to big scat¬ 
tering of the points in graphs caused by numerical inaccuracy. In this case we perform 
the averaging over all intersection points and proceed with this averaged point. Then we 
compute the average values of the inverse temperature 

i j<i 


and of the cumulant 


U. 


k 

av 


l 

To 




i j<i 
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FIG. 5: (Color online) The cumulant U\ for fc-th series as a function of /3 for 6 = 1 and various 
values of L x = 20, 40, 60, 80, 100: (a) Comparison of results for k = 1 (lines) and k = 2 (symbols). 
Note, that some pair of graphs (for example, for L x = 40 and L x = 80, 100 for k = 1) have no 
intersection within the considered range, (b) Results for k = 1 and intersection points (/?(), t/F). 


over ten intersection points for five system sizes for each set k. We observe (as expected), 
that average of crossing points do not coincide with crossing of average lines. 

Now we can compute the average values of the square of deviations from the mean values 
for the inverse temperature 


A/3 2 = 


1 


8 x 10 




k=1 


j<i 


and for the cumulant 


A U 2 a = 


8 x 10 




k=1 i j<i 

and evaluate numerical inaccuracy. In Fig. [6](a) we plot A/3 2 as a function of 6. The 
amplitude of the square of the inverse temperature deviation is about 10 -6 which is consistent 
with the average distance between points A/3 ~ 10 -3 . We do not observe pronounced 
minima in the interval 1.0 < 6 < 1.3, and the numerical inaccuracy is comparable with the 
scattering of the points. In Fig. [61(b) we plot data for the average square of the deviation for 
the cumulant AC/f. The point distribution is very similar to the ones obtained for the inverse 
temperature. The point for 9 = 1.1 is below its neighbors, but the distance between these 
points is of the order of the numerical inaccuracy. Therefore the above described procedure 
corroborates 9 — 1.1 as the optimal one. 
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FIG. 6: (Color online) (a) The average of the square deviation of the inverse temperature A/3 2 of 
crossing points as a function of 0; (b) The average of the square deviation of the cumulant AUf of 
crossing points as a function of 0. 

C. Computations of the susceptibility 


In this subsection we describe the results of the computation of the magnetic susceptibility 
at the critical point /3 C with various values of the exponent 9. Before we obtained two 
estimates for the critical /3 C : /3 ^ ax = 0.2565 (for method A) and /3,f ax = 0.2572 (for method 
B). In Figs. [71(a) and (b) we plot the magnetic susceptibility y as a function of the system 
size L x in the log-log scale for /3^ ax and (3 ^ ax , respectively. We expect the power law behavior 
y(L x ) oc LAj v for the more satisfying value of the exponent 9. In Fig. [7|(a) the deviation 
of the data from a straight line is less than in Fig. [71(b) therefore we select the estimate 
(3 C = /3^ ax = 0.2565 for future analysis. This is consistent with the analysis of Binder 
cumulants where it appears that /3,f ax seems to be out of the crossing window of U%. 

We perform a fit of the magnetic susceptibility by a linear function of L x in the log- 
log scale: x — a L b x with a and b being the parameters of the fit. Then we estimate the 
deviation of the values of susceptibility measured numerically from those obtained via the 
fitting function with the help of a chi-square y 2 defined as 

i— 1 


x . 2 = 


af 


(26) 


where N — 8 is the number of values Xi calculated for i = 1,2, ...,1V values of Lj = 
10, 20, 30,40, 50, 60, 80,100 the system size L X} f(x ) is the fitting function, of is the appro- 
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FIG. 7: (Color online) Magnetic susceptibility y as a function of the system size L x for various 
values of 6 = 1, 1.1, 1.2, 1.3, 1.4 at the point: (a) /3^ ax = 0.2565, (b) /3^ ax = 0.2572. 


priate variance defined by the error bars. 

The same number of MC steps is used for a given value of the exponent 9. Therefore, the 
variance u, is minimal for small values of the system size L x and increases for larger values 
of L x . The total number of MC steps decreases with increasing of 9 (from 5 x 10 6 for 9 = 1 
to 5 x 10 5 for 9 = 1.45). 

The parameter y 2 of the fitting procedure characterizes the “quality” of the data with 
respect to the proposed functional dependence. In our case this parameter describes the 
deviation of points from the straight line in the log-log representation. We plot y 2 as a 
function of 9 in Fig. [HI The general tendency is the decrease of y 2 with an increase of 9 , 
because due to the smaller number of MC steps the variance is larger in this case. We 
observe a minimum in the region 9 ps 1.2 — 1.3. Unfortunately, for this procedure we cannot 
evaluate the variance (inaccuracy) dy 2 . Therefore we repeat the procedure in more regular 
way. We split 320 impurity realizations (for each value of L x for a fixed 9) into 10 series 
and perform the fit by the formula ln(y) ~ ln(afc) + b & ln(L x ) for every series obtaining some 
value Xk- Then we average these values, compute numerical inaccuracy and plot the results 
in Fig. [9](a) (black circles for y 2 and red triangles for the average value). We can see, that 
y 2 reaches the minimum for 9* ~ 1.25. In Fig. [91(b) we plot the resulting power b as a 
function of 9. Value of b at 9* gives an estimate for r y/u± ~ 1.85 ± 0.05, using the window 
1.1 < 9 < 1.3. We present here also estimate 7/Vj_ = 1.90 ± 0.08 obtained from data of 
Fig. [9] for the value 9 = 1.1 recognized as an optimal value in subsection IV Bl at cumulant 
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e 


FIG. 8: (Color online) The x 2 parameter of the deviation from the linear fit (Eq. (1261) 1 of the 
magnetic susceptibility as a function of 9. 




6 0 

FIG. 9: (Color online) (a) The y 2 parameter of the linear fit with variance (xi for every series 
black circles, average result red triangles) as a function of (b) The result of the power of the fit 
b (bk for every series black circles, average result green triangles) as a function of 0. 

analysis. 


VI. CONCLUSION 

In this paper we have studied by MC simulations the scaling behaviour of thermodynam¬ 
ical quantities in the vicinity of critical point for 3d Ising system with randomly distributed 
parallel linear extended defects, modeled as non-magnetic impurities collected into lines 
along spatial direction z. We considered combined algorithm using Wolff and Metropolis 
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methods. Our results are consistently interpreted using the theory of anisotropic finite size 
scaling. 

We have estimated the value of the anisotropy exponent 6 using three different methods, 
namely, from dependence of correlation length £j_ on the linear size of the system near the 
critical point, from a temperature dependence of the fourth-order Binder’s cumulant, from 
finite-size scaling of the susceptibility. The values estimated are in the range 1.1 < 6 < 1.3 
and corroborate RG predictions for the model under consideration. 

We have also measured the value of 7/V_i_ from anisotropic finite-size scaling for the 
susceptibility. The value reported here y/Vj_ — 1.85 is a little bit below the corresponding 
RG estimates r y/u± ~ 2.0, r y/is± — 1.98. Value y/z/j_ ~ 1.9 estimated for 6 = 1.1 is in better 
agreement with theoretical results. We applied procedure described in this subsection to 
compute magnetization. However it does not give satisfactory results in this case. This is 
because magnetization is vanishing quantity at the critical point and therefore it is very 
sensitive to proper determination of critical temperature. 

Finally, let us mention that in spite of the substantial computation effort reported in this 
paper, the numerical values of the critical exponents are not extremely accurate. This is 


due 

Ref. 


;o the difficulty of the numerical techniques to deal with anisotropic systems (see e.g. 


521 ). but from the MC simulations presented in this paper, we believe that we can safely 


conclude in favor of the anisotropic critical point, since anisotropic scaling is overall nicely 
confirmed. 
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